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Abstract. In this paper, we study a strongly correlated quantum system that has 
become amenable to experiment by the advent of ultracold bosonic atoms in optical 
lattices, a chain of two different bosonic constituents. Excitations in this system are 
first considered within the framework of bosonization and Luttingcr liquid theory which 
are applicable if the Luttinger liquid parameters are determined numerically. The 
occurrence of a bosonic counterpart of fermionic spin-charge separation is signalled 
by a characteristic two-peak structure in the spectral functions found by dynamical 
DMRG in good agreement with analytical predictions. Experimentally, single-particle 
excitations as probed by spectral functions are currently not accessible in cold atoms. 
We therefore consider the modifications needed for current experiments, namely the 
investigation of the real-time evolution of density perturbations instead of single 
particle excitations, a slight inequivalence between the two intraspecies interactions 
in actual experiments, and the presence of a confining trap potential. Using time- 
dependent DMRG we show that only quantitative modifications occur. With an eye 
to the simulation of strongly correlated quantum systems far from equilibrium we 
detect a strong dependence of the time-evolution of entanglement entropy on the initial 
perturbation, signalling limitations to current reasonings on entanglement growth in 
many-body systems. 
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1. Introduction 

One of the key proposals in the field of quantum computing, made by Feynman, is 
the idea to use one quantum system to simulate another one in order to circumvent 
the problem of the qualitatively different complexity of quantum systems and classical 
computers as the standard simulation tools of science. The advantage of Feynman's 
approach lies in the fact that the new simulating quantum system may have advantages 
in experimental control: preparation, manipulation of parameters and measurement. 

Over the last decade, this proposal has been filled with life by the progress made in 
the preparation of dilute ultracold atom gases, highlighted by the now almost routinely 
achieved preparation of Bose-Einstein condensates. The use of Feshbach resonances or 
optical lattices has given us unprecedented control over interaction and dimensionality 
in strongly interacting quantum many-body systems of unique purity. This has been 
put to use in the creation of strongly correlated quantum systems that have been in 
the focus of interest in condensed-matter physics for a long timejl]. Examples are the 
observation of the superfluid to Mott insulator transition for Bose gases [2] and the 
fermionization of strongly interacting one dimensional bosons [Bl H]. 

But one can do more: it is possible to create physical systems of their own interest 
that have no counterpart in conventional condensed matter physics. To take one 
example, due to the internal spin degree of freedom of electrons it is quite natural in 
solids to consider models of two components of fermions. The existence of two fermionic 
components is the driving force of important phenomena such as collective magnetism. 

In this paper, we consider the bosonic equivalent of the two-component fermionic 
system. Strongly interacting two-component bosonic systems have no counterpart in 
condensed-matter physics, which is why they have found limited attention in the solid- 
state literature. Yet, they are quite easily implemented in the field of ultracold atom 
gases [HI E], and have generated quite some interest (for a review, see [7]). We focus 
on the case of one-dimensional two-component bosonic systems: on the one hand, they 
offer interesting physical phenomena such as a bosonic version [SI El HQ] of fermionic spin- 
charge separation which was discussed for cold atoms by [TTJ [TJl [131 [H] . On the other 
hand, very good control exists both analytically via bosonizationJTH] and numerically 
via (time-dependent) DMRGplDTj. 

In this paper, we start out by considering the low-energy physics which is 
characterized, as all other critical one-dimensional quantum systems, by a few effective 
Luttinger liquid parameters, which we determine both by a mapping in a limiting regime 
and more generally by DMRG. In principle, the Luttinger liquid parameters determine 
completely the static and linear response properties of the system. 

We move on to discuss the spectral functions which are the cleanest way to observe 
Luttinger liquid physics, in particular spin-charge separation and characterize the linear- 
response behaviour of two-component bosonic systems. These are obtained in the 
framework of dynamical DMRG. 

Experimentally, one will be confronted by certain limitations of ultracold atomic 
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systems: currently, spectral functions are unavailable and one is restricted to monitor 
the time-evolution of excitations. These in fact show a separation of symmetric and 
antisymmetric density combinations ("charge" and "spin") of the individual species 
from which results in good agreement with those from the spectral functions can be 
derived. These observation still hold if we also take into account that in current 
implementations of two-component bosonic systems there are slight differences in the 
intraspecies interactions of the two components. Moreover, the presence of a confining 
harmonic trap potential does not qualitatively alter the results, as we can show by 
explicit simulation using time-dependent DMRG. 

Last, but not least, we turn to the discussion of the time-evolution of the entropy of 
entanglement in the various out-of-equilibrium scenarios considered in this paper. This 
study was originally motivated by the fact that the efficiency of DMRG simulations is 
limited by entanglement growth which impacts exponentially on the numerical resources 
needed. Here, it turns out that there are surprisingly large variations in the generally 
accepted scenario of an essentially linear entanglement growth after quenches (as which 
all our scenarios can be interpreted). While we cannot give a general explanation for 
the phenomenon, we hope to provide a stimulus for further research to develop a more 
complete understanding of entanglement evolution. 



2. Model 



Cold atomic gases with two hyperfine species confined in optical lattices can be described 
by the two-component Bose-Hubbard model given by the following Hamiltonian[18j: 

H = —JJ2j,u \^j-\-i,v^3,v + h.c.^j + J2j,u J ' v 2 J V " /-q 

where j is a site index and v labels the two different flavours of bosons. J is the hopping 
strength, U u the intraspecies and Uu the interspecies onsite interaction. Here we assume 
that the two hyperfine species have the same mass and see the same lattice potential. 
Unless otherwise stated we use U\ — U 2 = U, and u = U/J, u 12 = U\ij J ■ £j, u describes 
an external potential, i.e. given by a trap. 

In the following we denote the lattice spacing by a. We mainly consider incommensurable 
fillings with equal densities n\ = n 2 = n smaller than one. For vanishing interspecies 
interaction we recover the one component Bose-Hubbard model being in a superfluid 
phase. The superfluid phase remains stable for finite U\2 up to U\2 ~ U. For 
U12 > U the interspecies interaction becomes dominant and a demixing of the flavours 
occurs, i.e. a phase separation[8j H9] . The physics corresponds then to the one 
of an itinerant ferromagnetic system with non-Luttinger liquid properties [20l [21]. 
Experimentally, to test for Luttinger liquid properties, the relevant regime is given 
by U\2 ~ U, with U12 slightly below U to avoid the phase-separation regime. Strongly 
differing interaction parameters have not been realized experimentally so far; our typical 
interaction parameters are chosen accordingly. 
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3. Approximations 

Continuum model In a weakly interacting superfluid phase in the low filling limit, the 
Bose-Hubbard model can be mapped to the continuous Lieb-Liniger model [2~2~j [23] . The 
mapping can be performed taking the limit a — > while leaving J a 2 constant. 
The Hamiltonian for two bosonic species in the continuum is 

+ V(x)*t(x)* u (x) + |(*t(x)) 2 (* v (x)) 2 ) 

+ ^fjdx (v&t(x)vl/ 1 (x))(vl/ t 2 (x)vl/ 2 (x)). (2) 

Here \E , ^- ) is the bosonic annihiliation (creation) operator, V the external potential, 
M the mass of the particles, and g and g^ are the strengths of the intra- and inter- 
species interaction, respectively. The parameters of the continuum model and the lattice 
model are related by J a 2 = ^7, the interaction strength to the 5- interaction strength 
by Ua = g and Vyia = gu, and the density p to the filling factor n by pa = n. In the 
hydrodynamic approximation [21] the sound velocities of this model are given by 

v c ,s = v / \]l±g 12 /g (3) 

with v = Jgp/M. (4) 

The indices stand for charge and spin, corresponding to symmetric (charge) and 
antisymmetric (spin) combinations of the two bosonic fields. At the symmetric point, 
where the interspecies interaction is equal to the intraspecies interaction, the model is 
SU(2) symmetric and can be solved by the Bethe ansatz. For the special case g\2 = g the 
sound velocity has been determined using the Bethe ansatz [211 [25]. The spin dispersion 
becomes quadratic at this point. 



Bosonization The low energy physics of the Bose-Hubbard model can be described 
by the bosonization approach[15]. Two bosonic fields, 6 U and <p u , are introduced 
in the continuum which are related to the phase and the amplitude of the original 
bosonic operator, respectively. To be more precise, in this representation the bosonic 
creation operator becomes b\,(x) = ^p^Y /p e l2p( - 7rpQX ~^^e~ ie '' and the density operator 
p u (x) = po — 1 / ttV 4> u (x) + p J2p^oe l2p ^ 7TpoX ~ (l '" i ' x ^- Here p is the average density, and 
-d x (j) v and 9 U are conjugate operators. The advantage of this representation is that 
the low-energy properties of the system are described by a quadratic Hamiltonian. 
Introducing the "charge" and "spin" degrees of freedom by <p c = 1/\/2(</>i + <p 2 ) and 
4> s = l/v / 2(0i — 02) it separates into two different part and is given by 

H =H C + H S with (5) 



H c = — I d.r 

2tt 



v c K c (d9 c ) 2 + ^(d x <b c 



and 

(6) 
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Vs K s {d6 s y + ^r{d x <j) s f 



+ / dx cos ^- ( ? ) 

Here the parameters v v are the velocities and K v the Luttinger liquid parameters for the 
corresponding field. Due to this mapping the asymptotic physics of the Bose- Hubbard 
model is totally determined by the velocities and the Luttinger liquid parameters. In 
the case of intermediate interaction strengths, the relations between the microscopic 
parameters of the Bose-Hubbard model and the parameters of the bosonization approach 
are difficult to establish analytically. In the limit of small interaction strength and 
low filling they are given for the single species Bose-Hubbard by Jamr = VqK and 
Ua/2n = v /K. 

In the limit of weak interspecies coupling the parameters for charge and spin degrees 
of freedom can be related to the parameters of the system without interspecies coupling 
by the perturbative expressions 



irK t 



V C , S = V ^l±^ 

K c>s = K / x fT± 9 -^ 







(8) 

Outside the range of validity of these two approximations the parameters can be 
determined numerically which we do in the following for a broad range of parameters. 
The Luttinger parameter is determined by using 

7T 7T 

K c = - v c k c , K s = - v s k s . (9) 

The compressibility [15J is given by k CjS = (^f^^J > where E (L) is the energy of the 
ground state for a system of length L and the combined "charge" and "spin" densities 
are denoted n c s — n\ ± n-i. The derivative can be computed numerically, yielding 

,.-1 ^ j E(N+AN,N+AN)+E(N-AN,N-AN)-2E(N,N) 

-1 ^ t E(N+AN,N-AN)+E(N-AN,N+AN)-2E(N,N) \ LV ) 
K s ~ -k 4AAf 2 

E(Ni, N 2 ) is the ground state energy of the system with Ni particle of the first flavour 
and N 2 of the second one. After a L — > oo extrapolation, we can obtain the Luttinger 
parameter by using Eq. (Q and the numerically determined velocities. Note that the 
formulae given for K s and k s in Equations and fllOp differ by (mutually cancelling) 
factors of 4 from the usual relationships in electronic models, where n s = (n\ — n 2 )/2; 
the factor 1/2 is not meaningful in the model considered here and hence dropped. The 
final value of K s is of course unaffected. 

Fig. [T] shows the dependence of the Luttinger parameters on the interspecies inter- 
acting strength u 12 . The charge Luttinger parameter K c decreases with increasing u± 2 , 
the spin Luttinger parameter increases. As expected, the bosonization result agrees with 
the numerical results for small u± 2 and starts to deviate with increasing ui 2 . The spin 
Luttinger parameter shows stronger deviations. Fig. [2] shows the Luttinger parameter 
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Figure 1. (Color online) Dependence of the charge and spin Luttinger parameter on 
the interparticle interaction strength u\2- A comparison of analytical results (line, see 
text) and numerical DMRG results (symbol) is shown. The parameters used are (a) 
u = 3, n w 0.63 and (b) u = 2, n w 0.88. 



for fixed interaction strength (u = 3, u\2 = 1.2, 2.7) and varying density n. Note that 
for the spin Luttinger parameter u^vq/ {txKq) is of the order of one, and thus the pertur- 
bative result (Eq. [8]) is not reliable anymore. The DMRG results deviates significantly 
from the perturbative results; the deviations increase with density. A similar effect also 
occurs in the one-component Bose-Hubbard model where for densities above 0.5 lattice 
effects becomes relevant, see the Appendix for more details. 



4. Spectral Functions 

The key quantity to characterize single particle excitations in many-body systems is the 
dynamic single-particle spectral function, because it can be probed easily in solid-state 
setups. The single particle spectral function is defined as 

A{q } LJ + i V ) = --Q{0\b q>1 — —. ^JJO), (11) 

7T Eq + UJ + IT] - H 



Excitations in two- component Bose-gases 



7 



u 12 = 1.2, charge 

• u 12 = 1 -2, spin « 
t u 12 = 2.7, charge 

♦ u 12 = 2.7, spin 
- - - Bosonization 



0.25 0.5 0.75 1 

n 

Figure 2. (Color online) Dependence of the charge and spin Luttinger parameter 
on the charge background density. A comparison of analytical results (line, see text) 
and numerical DMRG results (symbol) is shown. The parameters used are (a) u = 3, 
U12 = 1.2 and (b) u = 3, u\2 = 2.7. 

where |0) is the ground state with energy E . For fermionic two-component Hubbard 
model this function signals the separation of spin and charge degrees of freedom by two 
distinct peaks at different frequencies. 

In the following we determine this function for the two-species bosonic model using 
a variant of dynamical DMRG [26} [27] . Our formulation of the algorithm is entirely 
based on matrix product states. This allows to avoid the targetting of multiple states 
in the DMRG algorithm, which is achieved only at substantial numerical cost. The 
new formulation of the algorithm saves more than an order of magnitude in time. As 
always, DMRG prefers open boundary conditions leading to spurious effects in the 
spectral functions; in this context, filtering procedures have been used[26]. In order to 
reduce boundary effects we use the quasi-momentum definition of the Fourier transform 
b q ,v — J2j srn (f+i) &?> with the momentum \k\ = qir/(L + 1). These states create a 
single particle basis of eigenstates of the non-interacting system (u = U\i = 0) with 
open boundary conditions. They have nodes at the edges of the systems and are there- 
fore better suited to open boundary conditions than the usual definition of the Fourier 
transform using the eigenfunctions of the non-interacting system with periodic bound- 
ary conditions. 

Dynamical DMRG uses a correction vector method to calculate the spectral 
function: First define \lv(q,u)) = &ti|0). Then the correction vector is given by 
\cu(q,u + irj)) = w )) and thus obey 

(Eq + u + irj - H) \cv(q, to + in)) = \lv(q, u)). (12) 

This complex equation for the correction vector can then be solved using 
the GMRES algorithm. The value of the spectral function is thus A(q,u) = 
— -Q(lv(q,uj)\cv(q,uj- s riT])). A significant improvement is gained by considering the 
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derivative of A(q, uj + irj), 




(13) 



where the last line is only valid if \lv(q,uj)) = \lv(q,u)) and l^) denotes the complex 
conjugate of Using the derivative, we can determine the spectral function with a 
smaller number of correction vectors and thus much more efficiently. 

Note that in the following only the normalized spectral function 



is used, such that / du A norm (q, uj + irj) = 1 holds for every q. The full spectral function 
A n orm{q, is shown in Fig. [3j For the used system parameters L = 64 and r\ = 0.1 
we needed up to 2000 states for each correction vector. One clearly sees two different 
branches with a linear dispersion relation uj ~ v CjS q yielding two different velocities. 
Thus this system exhibits spin-charge separation. Fig. H] gives a more detailed view of 
the spectral function. For a number of given momenta q = kn/(L + 1) the spectral 
function is plotted versus ujq~ x . Therefore the norm of the scaled spectral function 
obeys / d (ujq^ 1 ) A norm (q, u + irj) = q~ l and the norm decreases with increasing q. The 
position of the peak is roughly given by uj ~ v C)S q, thus one expects two peaks at the 
two velocities v c and v s [28J . The spectral function shows a number of boundary effects 
similar to the one encountered for a single-component Bose-Hubbard model, see the 
Appendix for a detailed discussion of the one-component model. Here an addition peak 
at uj ~ occurs, visible for instance for q = 157r/65. Furthermore, for small values of q 
the charge peak is shifted to lower frequencies and above the charge peak an additional 
shoulder appears (i.e. for q = 107r/65). There are also some effects which do not occur 
for a single-component model: for large momenta the charge peak splits up and an 
additional peak below the main peak emerges. The plot also shows some effect of a non- 
linear dispersion. Both spin and charge peak are shifted to higher frequencies for higher 
momenta, similar to the effects visible for a Bose-Hubbard model, see Fig. IA21 In Fig. [3] 
momenta lower than q « 0.25 are not plotted, because for the system sizes considered, 
finite size effects dominate there and obscure the physically relevant information. 

5. Density perturbations and single-particle excitations 

In ultracold atom experiments, spectral functions are hard to observe, but using suitably 
tuned and focused lasers it is easy to create local density perturbations. Theoretically 
the density perturbations can be created applying an external potential of the form 
e j,cs ~ exp(— (J — jo) 2 I la j) for times t < 0. For times t > the potential is switched 




(14) 



off. 
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Figure 3. (Color Online) Density plot of the one-particle spectral functions A(q, u). 
The following parameters were used n = 0.625, u = 3, u\2 = 2.1 on a system with 
L = 64 sites and a broadening r\ = 0.1. The charge branch is above the spin branch. 
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Figure 4. (Color online) One-particle spectral functions A(q, ui) at momenta q = 
k/65ir/a plotted against ujq^ 1 . Two peaks corresponding to the spin and the charge 
excitation can be distinguished. The vertical lines mark the position of it CjS . The 
following parameters were used n — 0.625, u = 3, u±2 = 2.1 on a system with L = 64 
sites and a broadening rj = 0.1. 



Excitations in two- component Bose-gases 



10 




Figure 5. (Color online) Snapshots of the time-evolution of the charge and spin 
density distribution of a small charge (a) and spin density perturbation (b) created at 
time t = Oh/ J. The system parameters are 71^2 = 0.625, u = 3, U12 = 2.1. 



Time-evolutions can be calculated by adaptive time-dependent DMRG most easily 
by using Trotter decompositions of short-ranged Hamiltonians leading to local time- 
evolutions [291 EH EH]- For longer-ranged interactions or systems with large local state 
spaces, it is more efficient to consider global time-evolutions|32j. As we are dealing lo- 
cally with products of two bosonic state spaces, we calculated the time-evolution of the 
density perturbations numerically in the latter framework. The algorithm was formu- 
lated using matrix product states and the global time evolution was done using a Krylov 

algorithm for exponentiation [33] with a fixed error bound per timestep [31]. The used 

2 

error bound \^f(t + At)) — exp(— iHAt)\^(t)) is of the order of 10~ 5 with a timestep 
of At « 0.1 — 0.2. Usually, 6 to 10 Krylov vector were used. For the Krylov vectors up 
to 3000 states were used in the case of the density perturbations. The time-evolutions of 
single particle excitations are much harder to perform, up to 7500 states were used there. 

Snapshots of the density evolution of these excitations for one parameter set are 
shown in Fig. [51 One sees that the charge and spin perturbation which is created at 
time t — splits up into two counter-propagating excitations. The velocity of the spin 
perturbation is much lower than the velocity of the charge perturbation. Even after 
separating into two perturbations the amplitude shows a decay. This decay is very slow 
for weak inter-species interaction (see Fig. [5]). However, if the inter-species interaction 
is approaching the value of the intra-species interaction a strong broadening of the spin 
perturbation can be seen (see Fig. [6]). Due to the very low velocity and the broadening 
the two peaks only separate at very long times. The strong broadening of the spin peak 
is expected since in the limit of equal interaction strength the system has a quadratic 
spin dispersion relation. Therefore the region of the linear regime of the dispersion for 
u%2 < u decreases if u\2 gets closer to u. The dips in front of the spin perturbation and 
behind the charge perturbations are due to the remaining interaction between the spin 
and the charge degrees of freedom and finite-size effects. 



The results for a density perturbation connect directly to a setup where the time- 
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Figure 6. (Color online) Snapshots of the time-evolution of the charge and spin 
density distribution of a small charge (a) and spin density perturbation (b) created at 
time t = Oh/ J. The system parameters are 71^2 = 0.625, u = 3, U12 = 2.9. 



evolution of densities is followed upon the creation of single-particle excitations in 
bosonic gases instead of a density perturbation. In a one-dimensional system these single 
particle excitations decay into a charge excitation and a spin excitation. In contrast in 
a higher- dimensional system single particle excitations have a finite life-time. Since in 
cold atomic gases one can realize systems of different dimensionality, this second setup, 
if realized, would give the possibility to directly confront the decay in a one- dimensional 
system with the finite life-time in a three dimensional system. Fig. [7] shows the time 
evolution of the density and bipartite entanglement entropy profiles for a creation and 
an annihilation of a particle at time t = 0. The initial excitation splits up into a right 
and a left moving part and one can clearly see the two different velocities for the charge 
and spin density. This eventually leads to the separation into a spin and a charge 
excitation. The entropy profiles exhibit a significant difference between the creation 
and the annihilation of a particle. After removing one particle the entropy stays almost 
constant (up to some small waves) with respect to time, the entropy strongly increases 
between the two counterprogating peaks after adding an addional particle. The reason 
for this behaviour is unknown; observations and conjectures will be discussed in more 
detail in Sec. H 

6. Experimental constraints: Effects of unequal intraspecies interaction 
strengths and confining trap potentials 

An experimental realisation of a two-component Bose-Hubbard model is given by the 
use of two different hyperfine states of 87 Rb, for instance \F = 2, trip = — 1) and 
\F = 1, mp = 1). The s-wave scattering lengths for these states are approximatly 
Gt2 = 91.28ob, di = 100.4ob, where as is the Bohr radius [35]. Hence the intra- 
species interactions are slightly distinct for the two bosonic species. By comparison, 
the interspecies scattering length an is of the same order of magnitude and can be 
tuned by a Feshbach resonance pH [35]. Treating unequal intraspecies interactions in 
the bosonization approach results in a coupling of the spin and the charge part of the 
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Figure 7. (Color online) Snapshots of the time-evolution of the charge and spin 
density distribution of a single particle excitation ((a) - (c) creation of a particle, (d)- 
(f) annihilation of a particle) created at time t = Oft,/ J; (a), (d) at time t = Oh/ J, (b), 
(e) at time t = 1.5%/ J and (c), (f) at time t = 2.5h/J. The system parameters are 
m, 2 = 0.625, Ux/J = U 2 /J = 3., U 12 /J = 2.1. The charge density is shifted by 1.25 
for better visibility. 



Hamiltonian. To decide if for the experimentally relevant parameters the coupling is 
already important we have calculated the time-evolution of a single particle excitation 
of a system with realistic parameters. The parameters have been determined using a 
lattice depth of V = 4:.3E R and a interspecies scattering length a 12 = 80ae, where 
Er is the recoil energy of the optical lattice. This yields the following Bose-Hubbard 
parameter, UJJ = 2.983, U 2 /J = 2.712, U 12 /J = 2.377. Fig. E shows the density 
profiles for different times. Even though a small coupling between the spin and the 
charge degree of freedom might be present, the single particle excitation splits up into 
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Figure 8. (Color online) Snapshots of the time-evolution of the charge and spin 
density distribution of a single particle excitation created at time t = Oh/ J; (a) at time 
t = Oh/ J, (b) at time t = X.5H/J and (c) at time t = 2.5H/J. The system parameters 
are m, 2 = 0.625, Ux/J = 2.983, U 2 /J = 2.712, [/ 12 /J = 2.377. The charge density is 
shifted by 1.25 for better visibility. The arrows in (c) mark the clear separation of the 
charge and the spin density waves [Reproduced from Ref. |10j]. Copyright American 
Physical Society 
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Figure 9. (Color online) Snapshots of the time-evolution of the charge and spin 
density distribution of a single particle excitation created at time t = Oh/ J in a trap; 
(a) at time t = Oh/ J, (b) at time t = 1.5H/J and (c) at time t = 2.5h/J. The system 
parameters are m t 2 — 0.625, u\ — 3, U12 = 2.1 and ej jV = 6 • 10~ 3 (j — jo)" ■ The charge 
density is shifted by 1.25 for better visibility. 



a charge and a spin excitation. Therefore the separation in spin and charge survives a 
small experimental mismatch in the interspecies interaction strength. 

A further complication of the ultracold atomic gases setup is the presence of an 
parabolic trapping potential. In the Bose-Hubbard model this trapping potential can 
be described by adding = V (j — j ) 2 . 

We calculated the time-evolution of a single particle excitation with a trapping 
potential, see Fig. [9j In contrast to the case without a trap the ground state density 
is not constant anymore. Therefore the velocity of the charge and the spin excitation 
now depend on the spatial position of the excitation. Still, two counter-propagating 
waves can be observed and the spin-charge separation is not qualitatively changed. The 
effect of a trapping potential has already been discussed in the context of spin-charge 
separation for two component fermionic systems [TlT [T2| [13] , where similar robustness 
of the results was found. 
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7. Entropy of entanglement 

The superposition of states that is characteristic of quantum mechanics implies for 
many-body systems the phenomenon of entanglement that is the key deviation of the 
quantum from the classical world and a key resource of quantum computing. In this 
Section, we focus on one measure of entanglement in bipartite systems, the entropy 
of entanglement, which is defined in the case of pure quantum states for an arbitrary 
bipartition of the system into a left part A and a right part B by cutting at bond i by 
forming the reduced density operators of A and B, 



The entropy of entanglement is then given by the von Neumann entropy of either 
reduced density operator pA or pb- 



The identity of both definitions follows from the well-known observation that the 
eigenspectra of reduced density operators of a bipartition are identical. 

Entropy of entanglement has been studied in numerous contexts. In the present 
work, the focus is on the time-evolution of the entropy of entanglement in an out-of- 
equilibrium setting. On the one hand, this question is of fundamental interest for the 
understanding of coherent out-of-equilibrium quantum dynamics. On the other hand, 
it turns out that the time-evolution of the entropy of entanglement is closely related to 
the performance of time-dependent DMRG and related methods: the number of states 
(matrix dimensions) needed in such simulations are in a roughly exponential relationship 
with the entropy of entanglement. This means that e.g. a linear growth of entanglement 
entropy in time is reflected in an exponential growth in matrix dimensions, yielding the 
key limitation [36] for the time-scales accessible for such algorithms. 

Building on the Lieb- Robinson theorem [3 7j it has been shown by Osborne [38] 
that for arbitrary short-ranged Hamiltonians in one dimension, entanglement growth 
is bounded linearly in time, S(t) < S(0) + ct, reflecting a finite speed of propagation in 
such Hamiltonians, implying a potentially exponentially growth of matrix dimensions 
in time. In fact, such linear growth of entanglement has been observed [391 HOI E] and 
been traced back to the fact that ouf of equilibrium quantum states show excitations 
propagating through the system, leading to a linearly expanding "light cone" where 
information is exchanged between bipartition parts A and B. 

In this paper, we are considering three quite different types of time-evolution: 
(i) time-evolution after an insertion of a single particle, (ii) time-evolution after the 
extraction of a single particle, (iii) time-evolution after removing an external potential 
that created a density perturbation. All processes should generate excitations moving 
at the same maximal finite speed of propagation; one would therefore naively expect 
that all cases lead to a qualitatively similar linear growth of entanglement entropy in 
time, although we cannot expect identical results: for example, there is no particle-hole 
symmetry relating (i) and (ii). 





(15) 




(16) 
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Fig. [7] shows the time-evolution of the entropy profiles after a single particle ex- 
citation by either insertion (left) or removal (right), i.e. cases (i) and (ii) where the 
entanglement entropy has been measured for all possible bipartitions. The drops at the 
ends are a natural consequence of the small dimension of either p^ or ps there. The 
difference between (i) and (ii) is striking. In the case (i) entanglement entropy increases 
roughly linearly in time, as expected. However, in case (ii) entanglement entropy stays 
almost constant. This implies that simulating case (i) is much harder numerically than 
case (ii). For a comparison, we show case (iii), a small density perturbation evolving in 
time: the entanglement entropy is almost entirely unchanged under the time-evolution, 
see Fig. EH 

To summarise these observations, we show the maximal entanglement entropy on a 
chain versus time for cases (i)-(iii) in Fig. [TT], with the ground state entanglement given 
as a reference. As already seen in Fig. [7J a enormous difference between the creation 
and the annihilation can be observed. 

We can exclude that these observations are due to limitations of the numerical 
method that obviously neglects some of the entanglement entropy due to its inherent 
truncations. However, these results are converged in the used matrix dimensions; we also 
have calculated the fidelities between the initial states and the result of time-evolutions 
up to time t and then back to 0, finding fidelities that are essentially 1. The entanglement 
entropies shown can therefore be considered exact. 

In the case of a weak density perturbation, we attribute the observed behaviour 
to the fact that the out-of-equilibrium wave function is essentially a weakly distorted 
ground state wave function that does not affect the entanglement structure of the ground 
state substantially. It is much less obvious to explain the difference between (i) and (ii) 
despite the absence of a particle-hole symmetry. If we consider Fig. in the case of 
particle creation the increase of entanglement entropy is limited to the regions into 
which spin and charge excitations have already propagated that link ever larger regions 
of space. This is as expected. In the case of particle annihilation we might suppose that 
due to the relatively low density entanglement is removed because the chain is to some 
extent cut by the removal of a particle. In particular, if one thinks of the initial state as a 
superposition of many different Fock states, the application of the annihilation operator 
eliminates all states in which no particle occupies the site. Thereby the entanglement 
in the system is reduced. Indeed, a small dip is observed. The perturbation in the 
entanglement entropy propagates again linearly with the charge and spin excitations, 
but quantitatively it stays essentially unchanged. A possible explanation is suggested by 
perturbation theory. For short times the time-evolution operator can be approximated 
by an expression of the structure 1 + \tJY J ip\b i ±i + h.c.) + it{U /2)ni{rii — 1), where we 
have suppressed the multi-component nature of the problem for illustrative purposes. 
While the time-evolution operator applied far from the perturbation conserves the state 
and thereby the entanglement entropy, we see that close to the perturbation changes 
occur. At the densities considered, the site where the annihilation operator is applied has 
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Figure 10. (Color online) Entropy for a density perturbation n 12 = 0.625, Ui/J = 
U 2 /J = 3., I/ia/J = 2.1. 
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Figure 11. (Color online) Maximum of the bond entropy versus time for single particle 
excitations (both b^\0) and 6|0)) and for a small density perturbation. The small wave 
occuring for the b'\0) curve are an effect of the lattice. Note the enormous difference 
between the two single particle excitations. The system parameters are ni,2 = 0.625, 
u = 3, U12 = 2.1. 



most likely occupation number or 1 for the respective species in the contributing Fock 
states. After application of the operator, the Fock states either vanish or have occupation 
number on the site, such that only one of the terms in the hopping operator couples 
to the state at the site of the perturbation. Rerunning the argument for the creation 
operator, no contribution vanishes (if we assume intermediate interaction strengths) and 
all hopping terms couple, which might indicate that changes in entanglement are much 
more pronounced in the latter observed. This is obviously not rigorous at all, and 

at the moment we can only conclude that the simple picture of excitations transporting 
entanglement leading to linear entanglement growth needs substantial refinement. 
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Appendix A. Single-species Bose-Hubbard model 

The single-species Bose-Hubbard model has been discussed extensively in previous 
literature; for an overview of the literature, see pQ. In this appendix we would like 
to demonstrate some technical details useful to understand our results for the two- 
species Bose-Hubbard model in more detail at the example of the one-species model. 
The Hamiltonian of the one-species Bose-Hubbard model is given by 

H= -J^(b) +l b h + h.c) + Y. 3 U -^^ (A1) 

+ J2j e j n ji 

where we used the same notation as for the two-component model. Let us first comment 
on the deviation we see comparing the values of the velocities and the compressibility 
with the analytical formulae for the continuum limit. For the single component Bose- 
Hubbard model the expressions are given by 

K Q = K Q / (ttv ) 

K = ^(i-^/(27r))"* (A.2) 

v = 2rV7(l-V7/(27r)) 1 

with interaction parameter 7 = U/2n. In Fig. IA1I we compare the numerical results for 
the compressibility Kq = Kq/{ttvq) with the analytical formula for the continuum limit. 
The Lieb-Liniger results yield a good approximation for the compressibility for small 
densities n. This corresponds to rather high values of 7 of the order of 10. In contrast, 
for larger densities above 0.6 the results deviate considerably even for small values of 
7 « 1. Kollath et al. showed in [42] that the velocity of small density perturbation of 
the Bose-Hubbard model is given by the Lieb-Liniger solution for n < 1 and 7 up to 4. 
Therefore we see the same effect as in the case for the two-component model, that the 
deviations for the compressibility at n ~ 1 become larger than the deviations for the 
velocities. 



To investigate finite-size effects on the single-particle spectral function we show 
here the results for the single-component Bose-Hubbard model. Two different system 
lengths were considered in order to explore the boundary effects of a finite-size system. 
In fig. IA2l we show the spectral function A(q, uS) as a function oiu/q. Assuming a linear 
dispersion relation, one would expect a peak of the spectral function at the velocity v , 
see |43j. The numerical results indeed show this behaviour. Furthermore, a number of 
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Figure Al. (Color online) Compressibility k = K/vtt of a one-component Bose - 
Hubbard model for u — 2 and u = 3. 
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Figure A2. (Color online) Finite size scaling of the one-particle spectral function 
A(q, lo) at momenta q — 5/657r/a and q — 10/1297r/a respectively for a one-component 
Bose-Hubbard model. The vertical lines mark the position of uq. The following 
parameters were used n = 0.625, u — 3 on a system with L — 64 and L — 128 
sites and a broadening 77 = 0.1. 



boundary effects are visible: At ~ there exists an additional peak and the main 
peak is shifted to lower frequencies. Above the main peak there is also an additional 
shoulder. These are the same effects as we find for the two-component model. Let us 
note that all boundary effects become less pronounced for larger values of q. 
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